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Abstract 

Understanding how the dynamics of a neural network is shaped by the network structure, 
and consequently how the network structure facilitates the functions implemented by the 
neural system, is at the core of using mathematical models to elucidate brain functions. This 
study investigates the tracking dynamics of continuous attractor neural networks (CANNs). 
Due to the translational invariance of neuronal recurrent interactions, CANNs can hold 
a continuous family of stationary states. They form a continuous manifold in which the 
neural system is neutrally stable. We systematically explore how this property facilitates the 
tracking performance of a CANN, which is believed to have clear correspondence with brain 
functions. By using the wave functions of the quantum harmonic oscillator as the basis, 
we demonstrate how the dynamics of a CANN is decomposed into different motion modes, 
corresponding to distortions in the amplitude, position, width or skewness of the network 
state. We then develop a perturbative approach that utilizes the dominating movement of 
the network's stationary states in the state space. This method allows us to approximate the 
network dynamics up to an arbitrary accuracy depending on the order of perturbation used. 
We quantify the distortions of a Gaussian bump during tracking, and study their effects 
on the tracking performance. Results are obtained on the maximum speed for a moving 
stimulus to be trackable and the reaction time for the network to catch up with an abrupt 
change in the stimulus. 
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1 Introduction 



The brain performs computation by updating its internal states according to dynamical rules 
that are determined by both the nature of external inputs and the structure of neural networks. 
A neural system acquires its network structure through either learning from experience, or 
inheriting from evolution or both. Its impact on network dynamics is twofold: on one hand, 
it determines the stationary states of the network that lead to associative memory; and on the 
other hand, it carves the landscape of the state space of the network as a whole, which may 
contribute to other cognitive functions, such as movement control, spatial navigation, population 
decoding and object categorization. Understanding how the dynamics of a neural network is 
shaped by its network structure is at the core of using mathematical models to elucidate brain 



functions (Dayan and Abbott 2001). 



Associative memories, such as those based on the Hopfield model (Hopfield 1984; Domany, van Hemmen, and S( 
have been intensively studied over the past decades. However, an equally important issue has not 
been investigated to the same extent, namely, how a large-scale structure of the state space of a 
neural system may facilitate brain functions. Recently, a type of recurrent network, called contin- 
uous attractor neural networks (CANNs), has received wide attention (jWilson and Cowan 1972} 



lAmari 1977HGeorgopoulos, Kalaska, Caminiti, and Massey 1982|IMaunsell and Van Essen 1983t 



Funahashi, Bruce, and Goldman- Rakic 1989 Wilson and McNaughton 1993 Rolls, Robertson, and Georges-Fran( 



Ben-Yishai, Lev Bar-Or, and Sompolinsky 1995 Zhang 1996 Seung 1996 Samsonovich and McNaughton 1997 



Camperi and Wang 1998tlErmentrout 1998HHansel and Sompohnsky 1998[ITaube 1998HDeneve, Latham, and Po 



|Wang 200T Laing and Chow 2001 Wu, Amari, and Nakahara 2002 Stringer, Trappenberg, Rolls, and Aranjo 20( 



Brody, Romo, and Kepecs 2003 Erlhagen 2003 Gutkin, Pinto, and Ermentrout 2003 Renart, Song, and Wang 2( 



Trappenberg 2003t IFolias and Bressloff 20041 IWu and Amari 20051 IChow and Coombes 20061 



IMiller 2006} Machens and Brody 2008 Wu, Hamaguchi, and Amari 2008). These models are 



particularly useful when modeling the encoding of continuous stimuli in neural systems. Different 
from other attractor models, CANNs have a distinctive feature that their neuronal interactions 
are translationally invariant. As a result, they can hold a family of stationary states that can be 
translated into each other without the need to overcome any barriers. In the continuum limit, 
these stationary states form a continuous manifolcj^] in which the system is neutrally stable, and 



manifold is a mathematical concept. Simply stated, it is a set of points with a coordinate system. Here, 
each stationary state of the network is a point in the manifold. As will be shown later, these stationary states 
have the shape of a Gaussian bump, and their peak positions define the coordinates. 
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the network state can translate easily when the external stimulus changes continuously. Beyond 
pure memory retrieval, this large-scale structure of the state space endows the neural system 
with a tracking capability. An illustration of the landscape of the state space of a CANN and 
how it leads to neutral stability of the network stationary states is shown in Fig. 1. 

The tracking dynamics of a CANN has been theoretically studied by several authors in the lit- 
erature (see, e.g.,(Ben-Yishai, Lev Bar-Or, and Sompolinsky 1995 Zhang 1996t|Samsonovich and McNaughton 19 
Xie, Hahnloser, and Seung 2002^ IFolias and Bressloff 200"1 IWu and Amari 20051) ^. These stud- 



ies have demonstrated clearly that a CANN has the capacity of tracking a moving stimulus 
continuously and that this tracking property can describe many brain functions well. Detailed 
rigorous analyses of the tracking behaviors of a CANN are still lacking, however, including, for 
instance, 1) the conditions under which a CANN can track a moving stimulus successfully; 2) 
how the network state is distorted during the tracking; and 3) how these distortions affect the 
tracking performance of a CANN. In this study, we systematically investigate these issues. 

We use a simple, theoretically solvable model of a CANN as the working example. The sta- 
tionary states of the network have a Gaussian shape, and the network dynamics can be described 
using the basis functions of quantum harmonic oscillators. These analytical results allow us to 
explore the tracking performances of a CANN in much more detail than in previous studies based 
on simulation observations. In particular, we demonstrate clearly how the dynamics of a CANN 
is decomposed into different modes corresponding to distortions in the height, position, width or 
skewness of the network bump states. Their contributions to the tracking performance are de- 
termined by the amplitudes of the corresponding eigenvalues of the neuronal interaction kernel. 
For a CANN, neutral stability implies that the mode of positional shift dominates the network 
dynamics and that other modes have high-order contributions. From this property, we develop 
a time-dependent perturbative approach to simplify the network dynamics. Geometrically, this 
corresponds to projecting the network dynamics on its dominating motion modes. The solution 
of the perturbation method is expressed in a simple closed form, and we can approximate the 
network dynamics up to an arbitrary accuracy depending on the order of perturbation used. 
With this method, we explore two tracking performances of a CANN, namely, the condition 
under which the network can successfully track a moving stimulus and the reaction time for the 
network to catch up with sudden changes in the stimulus. Both properties are associated with 
the unique dynamics of a CANN and can be tested in practice. 

The idea of utilizing the special structure of a CANN to simplify its dynamics has also been 
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reported in other studies. Wu et al projected the dynamics of a CANN on its single dominating 
motion mode, the tangent of the attractor space, and simphfy the network dynamics into a one- 



dimensional equation (Wu, Hamaguchi, and Amari 2008). Several other authors also proposed 



different approaches to simplifying the network dynamics (Ben-Yishai, Lev Bar-Or, and Sompolinsky 1995 



IWennekers 2002 1 Renart, Song, and Wang 2003 : IMiller 2006]) . Compared with these studies, our 
contributions are that: 1) we elucidate the tracking behavior of a CANN in much more detail, 
particularly the effect of distortions in the network state on the tracking performance and the 
development of a perturbative method to systematically improve the description of the network 
dynamics; 2) we predict on the tracking behavior of a CANN, namely, the maximum speed 



for a moving stimulus to be trackable by the network (Hansel and Sompolinsky 1998) and the 
logarithmic nature of the reaction time; and 3) we also analytically investigate the tracking 
dynamics of a two-dimensional CANN. We expect that the mathematical framework developed 
in this work will have wide applications in the theoretical study of CANNs. 

The organization of the paper is as follows. In Section 2, the unique dynamics of a CANN is 
explored, and we demonstrate how a CANN is decomposed into different modes. In Section 3, 
a simple method of projecting the network dynamics onto its single primary motion mode, i.e., 
the tangent of the manifold formed by all stationary states, is used to investigate the tracking 
performance of a CANN. In Section 4, a perturbative approach is developed to systematically 
improve the description of the tracking dynamics of the network. The effects of the shape 
distortions of the network state on the tracking dynamics are described. In Section 5, the 
dynamics of a two-dimensional CANN is investigated. In Section 6, we present the discussion 
and conclusions. The mathematical details are presented in the appendices. Preliminary results 



have been pubhshed in ( Fung, Wong, and Wu 2008 ) . 



2 Intrinsic Dynamics of a CANN 

Let us consider a one-dimensional continuous stimulus x encoded by an ensemble of neurons. 
For example, the stimulus may represent the direction of movement, the orientation or a general 
continuous feature of objects extracted by the neural system. We consider the range of possible 
values of the stimulus being much larger than the range of neuronal interactions. We can thus 
effectively take x G (— oo,cxd). 

Let U{x, t) be the synaptic input at time t to the neurons whose preferred stimulus is x, 
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Figure 1: The energy landscape of the state space of a CANN. A canyon is formed by projecting 
the stationary states onto the subspace formed by 6Ui and 6Uq, i.e., the position and the height 
distortions of the network states. Motion along the canyon corresponds to the positional shift 
of the bump (inset). 



and let r{x,t) be the firing rate of these neurons. r{x,t) increases with the synaptic input, but 
saturates in the presence of global activity-dependent inhibition. A solvable model that cap- 



tures these features is given by the divisive normalization (Deneve, Latham, and Pouget 1999 

C/(x,t)2 



Wu, Amari, and Nakahara 2002] ): 

r{x, t) 



1 



l + kpJ:^dx'U{x',t)^' 
where p is the neural density and is a small positive constant controlling the strength of global 
inhibition. 

The dynamics of the synaptic input U{x,t) is determined by the external input /cxt(a;,i), 
the network input from other neurons, and its own relaxation. It is given by 

dU{x,t) 



dt 



/oo 
dx'J{x, x')r{x' ,t) — U{x, t), 
-oo 



(2) 



where r is a time constant, which is typically on the order of 1 ms (jErmentrout 199"8t Gutkin, Pinto, and Ermentrou 



Wu, Hamaguchi, and Amari 2008), and J{x,x') is the neural interaction from x' to x. 

The key character of CANNs is the translational invariance of their neural interactions. In 
our solvable model, we choose Gaussian interactions with a range a, namely. 



J{x, x') 



2'Ka 



exp 



t\2 



{x — x') 
2^ 



(3) 



where ^ is a constant. 

CANNs with other neural interactions and inhibition mechanisms have been studied in the ht- 
erature (see, e.g., (lAmari 1977| Ben-Yishai, Lev Bar-Or, and Sompohnsky 1995 lErmentrout 1998t 
Wang 200ip ). Here, we choose this model for the advantage that it permits an analytical solu- 



tion of the network dynamics. Nevertheless, the final conclusions of our model are qualitatively 
applicable to general cases. 

We first consider the intrinsic dynamics of the CANN model in the absence of external 
stimuli. It is straightforward to check that, for < k < kc = yl^/?/(8\/27ra), the network holds 
a continuous family of stationary states, which are (see Fig. 2A) 



U{x\z) 
f{x\z) 



Uq exp 



ro exp 



4a2 

(x - z)^ 
2^ 



(4) 
(5) 



where C/q = [1 + (1 - k/kcy/^]A/{4:^ak) and ro = [1 + (1 - k/kc)^/'^]/{2V27rakp). These 
stationary states are translationally invariant among themselves and have a Gaussian shape 
with a free parameter z indicating their positions, z is therefore the coordinate of the manifold 
formed by all stationary states. 

The stability of the Gaussian bumps can be studied by considering the dynamics of fluctu- 
ations. Consider the network state U{x,t) = U{x\z) + 6U{x,t). As derived in Appendix B, we 
have 

8 f°° 

T—5U{x,t)= dx'F{x,x'\z)6U{x',t)-5U{x,t), (6) 



where the interaction kernel F{x,x'\z) is given by 

2pU{x') 



F{x, x'\z) 



B 



J{x,x') — kp / dx" J{x,x")r{x") 



(7) 



We are interested in the eigenfunctions and eigenvalues of the kernel F{x,x'\z). Since the 
kernel is not symmetric with respect to the permutation of x and x' , we have to distinguish 
between its left and right eigenfunctions. To compute them, we need to choose a convenient 
set of basis functions. They should satisfy the boundary condition that they vanish when x 
approaches ±00. A convenient choice is one that contains a basis function taking the shape of 
a Gaussian bump. Since the quantum ground state of the harmonic oscillator potential has a 
Gaussian shape, a natural choice of basis functions is the wave functions of quantum harmonic 
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oscillators (jGriffiths 2004p . namely, 



Vn{x\z) 



1 



: exp 



(x — z)^ 
4a2 



X — z 



V2o 



(8) 



^(27r)i/2an!2 

where are the nth-order Hermite polynomials for n = 0, 1, 2, . . ., which can be expressed as 
the derivatives of the Gaussian function (|Griffiths 2004P 



, , , (-l)"(^/2a)"-V2 

Vn{x Z) = / , exp 



(x — z)^ 
4a2 



dx 



exp 



(X — 



(9) 



This choice of basis functions has the further advantage that an abundance of mathematical 
properties are available for our analysis, as summarized in Appendix A. 
Fig. 2B pictures the first four wave functions. We see that: 



The ground state of the quantum harmonic oscillator, 

U{x\z) 



vo{x\z) 



C/ov/(2vr)V2a 
corresponds to the stationary state of the network. 

The first excited state of the quantum harmonic oscillator, 

2a dU{x\z) 



(10) 



fill 



vi{x\z) = - 

Uo^{27r)y^a ^z 

corresponds to the positional shift of the stationary state. 

The second and third excited states, V2{x\z) and v^{x\z), respectively correspond to the 
distortions in the width and skewness of the bump. 



To express 5U{x, t) in terms of the eigenfunctions Vn{x\z) using Eq. (6), we consider F(x, x'\z) 
as an operator acting on Vnix\z). Since Vnix\z) is a basis, we have 

/oo 
dx'F{x,x'\z)Vnix'\z) ='^nmix\z)Fmn, (12) 

where F^n is the representation of the kernel -F(x,x'|2:) in the basis Vn{x\z), given by 

/oo rco 
dx / dx'vm{x\z)F{x,x'\z)vn{x'\z). (13) 
-oo J ~oo 
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Figure 2: (A) The Gaussian-shaped stationary states of the network. (B) The first four basis 
functions of the quantum harmonic osciUators, representing distortions in the height, position, 
width and skewness of the Gaussian bump. 



As derived in Appendix B for the particular choice of the output function ([T]) and the neural 
interactions ([3]), the elements of the matrix Fmn are given by 



TP 



1 - v/1 - k/kc, 

)l-n / n\ (-1)" 



71 — m / ^ 1 , 

I I ' 



0, 



m = n = 0; 
, n — m being an even integer; 

otherwise. 



(14) 



By using the matrix -Fmn, we can calculate the eigenvalues and the first four right eigenfunc- 
tions of the kernel F{x,x'\z), which are (see Appendix B) 



Ao = 1 - \Jl- k/kc, 



(15) 



An 

uf{x\z) 
u^{x\z) 

uf{x\z) 



2^-", for n> 1,..., 

vo{x\z), 

vi{x\z), 
n/T72 



^0 

1 



l-2Vl-fc/fc, 
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where 



^vi{x,z) + y -^3(2;, 2:), 



/ V V^' 

l-2^l-k/kA +1/2 



(16) 
(17) 
(18) 
(19) 

(20) 
(21) 



We note that Ui{x\z) = vi{x\z) and Ai = 1. Fig. 3A illustrates the first four right eigenfunctions 
of F. 

The right eigenfunctions of F correspond to the various distortion modes of the bump. 
For example, Uq{x\z) corresponds to the amplitude distortion of the bump, and Ui{x\z) to 
the positional shift. Starting from n = 3, u^{x\z) are linear combinations of Vk{x\z), where 
k = n, n — 2, . . .. Some of them may be dominated by a specific type of distortion, such as 
in Eq. (j2Up being dominated by the skewness. However, generally speaking, various distortion 
features have been coupled in the eigenfunctions of the interaction kernel. To avoid confusion, 
we continue to use Vn{x\z) as our basis functions in the rest of the paper. 

To identify the contributions of u^{x\z) to the network dynamics, we express 6U{x,t) in 
([6]) as a linear combination of u^{x\z), namely, 5U{x,t) = J2n^^n{z,t)u^{x\z). Using the 
orthonormality of the left and right eigenfunctions of F{x,x'\z), we have 



dx5U (x, t)u^ ix\z). 



(22) 



Assume that the motion of the bump is slow, so that dz/dt becomes negligible in Eq. ([6]); as we 
shall see, this assumption is valid as long as the external stimulus is sufficiently weak. Then, 
the projection of Eq. ([6|) on the eigenfunctions become 



T^6Un{z,t) = (A„ - l)6Un{z,t). 



Hence, 



6Uniz, t) = 6Uniz, 0) exp 
where SUn{z,0) is the initial value of the projection. 



(1 - Xn)t 



(23) 



(24) 



Eq. tells us that, since Ai = 1, the bump distortion 5Ui, which corresponds to the 
positional shift of the bump, is sustained under the network dynamics. This is the mathematical 
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expression of the neutral stability of the stationary states of the network. The eigenvalues for 
all other eigenfunctions arc all less than one. Since these other modes of bump distortions decay 
exponentially with the time constant t/(1 — An), their contributions on the network dynamics 
are small when compared with that due to the positional shift. An intuitive illustration of the 
decomposition of the network dynamics is given in Fig. 3B. 




+ 



Figure 3: (A) The first four right eigenfunctions of the interaction kernel F{x,x'\z) with 
k/kc = 1/2. (B) The distortions of the network state are decomposed into three parts, namely, 
distortions in the height, position and all other shapes. 

The fact that the dynamics of a CANN is dominated by only a few motion modes suggests 
that we may use this property to simplify the description of the network dynamics; that is, 
we can project the network dynamics onto its dominating modes. By this, we can analytically 
quantify the tracking bahviors of a CANN. 
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The Landscape of the State Space 

It is instructive to illustrate the energy landscape in the state space of a CANN. In the dynamical 
system that we are studying, we do not have a Lyapunov function. This is evident from the 
fact that the interaction matrix Fmn is not symmetric. Nevertheless, at each peak position z, 
one can define an effective energy function E{z) = J^ni^ ~ ^n)[SUn{z)]'^ /2, with 6Un{z) being 
the projection of U{x) — U{x\z) on the eigenfunction Un{x\z). Then, the dynamics in Eq. ([6]) 
or (I23p can be locally described by the gradient descent of E(z) in the space of 6Un{z). E[z) 
has the minimum value of zero when 6Un{z) = for n ^ 1 and 5Ui{z) can take any value since 
Ai = 1. This corresponds to the picture that the bump's position shifts, but the shape of the 
bump remains unchanged. When z varies, the solution of E{z) = traces a curve in the state 
space, facilitating the local gradient descent dynamics of the network, with the tangent of the 
curve at each point defining the direction of neutral stability. We can envisage this by a canyon 
surrounding the curve as shown in Fig. 1, and a small force along the tangent of the canyon can 
move the network state easily. 

The External Input 

To proceed, we need to define the external input. Without loss of generality, we choose 

(x - ZoY 



+ ar]{x,t), (25) 



4a2 

where a and a indicate the strength of the signal and the noise, respectively. r]{x,t) is Gaussian 
white noise, satisfying {rj(x,t)) = and {r]{x,t)ri{x' ,t')) = 6{x — x')5{t — t'). The signal term 
has the effect of driving the bump to the stimulus location zq, and the noise term tends to shift 
the bump position randomly. Our results are qualitatively independent of the exact form of the 
external input. We may choose other forms for the external input. The main conclusions will 
be the same once the chosen external input has a unimodal shape and has the effect of driving 
the bump to the stimulus position. 

The Setting of the Simulation Experiments 

We will use simulation experiments to confirm the theoretical analysis. Here, we describe the 
setting of simulation experiments to be used below. The network consists of N neurons. Their 
preferred stimuli x are evenly distributed in the range (— tt, vr] and satisfy the periodic condition, 
namely, U{x,t) = U{x + 27r,t). In the simulation, the integration in the continuum limit is 
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changed to the summation accordingly, i.e., p J dx —>■ J2i=i- We choose a = 0.5 <C 2n, so that 
comparison with the theoretical analysis for x £ (—00,00) holds. 



3 Tracking Dynamics in the Weak Input Limit 



We consider first the case when the external input is sufficiently small, i.e., a <^ 1 and o" ^ 1. In 
this case, we can reasonably assume that the network dynamics is dominated by the positional 
shift of the bump, and we can thus ignore other motion modes, that is, U{x,t) ~ U{x\z{t)), 
with only the bump position z{t) varying with time. As will be shown in Section 4, this actually 
corresponds to n = perturbation in the general perturbative method. 

After projecting the network dynamics, Eq. ([2]), on the basis functions vi{x\z{t)), we obtain 
(see Appendix C) 



dz 
lit 



a 



[z - zojexp 



8a2 



+ V2D7]iit), 



(26) 



where D = \f2aa'^ j {\JqT^ ^) is the diffusion constant, and r]\(t) is the overlap of r]{x,t) with 
vi{x\z{t)) given by 

POO 

ryi(t) = / dxr]{x,t)vi{x\z{t)), (27) 



which satisfies 



(r/i(t))=0 and (r/i(t)?7i(t')) = 5(t - *')• 



(28) 



This equation represents a one-dimensional Ornstein-Uhlenbeck process for the motion of 
the bump position in the state space ()Tuckwell 1988^ . Its meaning is straightforward: the 
first term on the right-hand side of the equation represents the contribution of the stimu- 
lus, whose effect is to pull the bump to the stimulus position zq; the second term represents 
noise, which tends to shift the bump position randomly. This equation agrees with Eq. (3.7) in 
(Wu, Hamaguchi, and Amari 2008), which considers \z — zqI <C a and hence ^-(^-^0)^/80.^ ~ i 



The distribution of the bump position can be calculated. The pull term in Eq. ([26]) can be 
written as the gradient of a potential. 



V{z) 



■ exp 



{z-zoY 



8a2 



(29) 



so that Eq. ()26p can be written as a Langevin equation, 

dz dV r-F^ , , 



(30) 
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This implies that z satisfies a Boltzmann distribution ()Ma 1976p . 




(31) 



Since the pull is maximum at z = zq i: 2a, two noise regimes can be identified. In the low 
noise regime, where {{z — zqY) <^ Aa?, the bump position is effectively trapped in a parabolic 
potential, whereas in the high noise regime with ((z — zq)^) S> 4a^, the pull vanishes and the 
motion of the bump resembles a random walk with {[z{t) — zqY) = 2Dt loosely coupled to the 
stimulus. 

However, simulation results show that Eq. (I3ip underestimates the variance of the bump 
position. This is because we have considered only the overlap of ri(x,t) with vi{x\z{t)). As will 
be shown in the perturbative approach, overlaps with other components cannot be neglected. 

By using the simplified dynamics in Eq. ()26|) . we explore two tracking behavior of the network, 
namely, the tracking of a moving stimulus and the catching up of an abrupt change in the 
stimulus. 

3.1 The Condition for Successful Tracking 

Suppose that the stimulus is moving at a constant velocity v and the noise is low and negligible. 
The dynamical equation becomes identical to Eq. (I26p after the transients, with zq replaced by 
vt. Denoting the lag of the bump behind the stimulus as s = zq — z, we have 



The size of the lag s is determined by two competing factors: the first term represents the 
movement of the stimulus, which tends to enlarge the lag; and the second term represents 
the collective effects of the neuronal recurrent interactions, which tends to shorten the lag. The 
tracking is achieved when these two forces match each other, i.e., v = g{s); otherwise, s diverges. 

Fig. 4A shows that the function g{s) is concave and has the maximum value of ^max = 
2aa/{Ty/e) at s = 2a. This means that if f > gmax; the network is unable to track the stimulus. 
Thus, (7max defines the maximum trackable speed of a moving stimulus to the network. Notably, 
^max increases with the strength of the external signal a and the range of neuronal recurrent 
interactions a. This is reasonable since it is the neuronal interactions that induce the movement 



ds 
dt 



dzQ dz 
dt dt 



as 

= V exp 

r 





(32) 
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of the bump, gmax decreases with the thne constant of the network r, as this reflects the 
responsiveness of the network to external inputs. 

On the other hand, for v < gma.x, there is a stable and unstable fixed point, denoted as 
si and S2 respectively, in Fig. 4A. When the initial lag is less than S2, it will converge to si. 
Otherwise, the tracking of the stimulus will be lost. Intuitively, this means that the propagation 
of neuronal activities in the network needs to be quick enough to compensate for the initial lag 
behind the stimulus. Figs. 4B and 4C show that the analytical results of Eq. (|32p agree well 
with the simulation results. 

3.2 The Reaction Time 

Another tracking behavior we study is the response to a stimulus experiencing an abrupt change. 
We measure how long it takes for the network to catch up with this change, i.e., the reaction time. 
Suppose that the network has reached a steady state at t < 0, and the stimulus position jumps 
from to zq suddenly at t = 0. This is the typical setting in experiments studying mental rotation 



behavior (Georgopoulos, Taira, and Lukashin 1993; Ben-Yishai, Lev Bar-Or, and Sompolinsky 1995 



IWu and Amari 2005p . The reaction time can be obtained by integrating Eq. (|26p with the initial 
value z{t = 0) = 0. 

When the noise is low and the jump size \zo\ is much smaller than the range a of neuronal 
interactions, we can derive an analytical expression for the reaction time. In this case, Eq. (|26p 
is approximated as 

dz 

T— « a{zo - z). (33) 
at 

We define the reaction time T as the time it takes for the bump to move to a small threshold 
distance 6 from the stimulus. T is calculated to be 

r.lM(M). ,34) 

This reveals that with small jump sizes, the reaction time of a CANN increases logarithmically 
with the jump size. 

Fig. 5 compares the simulation results with our theoretical predictions based on Eq. (j26p . 
We see that with small jump sizes, they agree very well with each other, and the reaction time 
indeed has a logarithmic nature in its increase with the jump size. However, we also notice 
that with large jump sizes, there is considerable discrepancy between the simulation results and 
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the theoretical predictions. This is because we have only considered the positional shift in the 
network dynamics. To improve the results, more motion modes must be included. 



4 Tracking Dynamics: A Comprehensive Study 

In this section, we develop a time-dependent perturbation approach to systematically improve 
the description of the network dynamics. Since the network dynamics is primarily dominated 
by the positional shift of the bump, with high-order distortions from other motion modes, we 
consider the network states with the following form, 

oo 

U{x, t) = U{x\z{t)) + an{t)vn{x\z{t)), (35) 

n=0 

where the bump position z{t) and the coefficients {a„(t)} specify the network state. 
The external input can be similarly expressed as 

oo 

4xt(a;,t) = Init)vn{x\z), (36) 

n=0 

with /„ = J dxIc:>,tix,t)Vnix\z). 

By substituting Eqs. (j35p and (|36p into Eq. and making use of the orthonormality and 
completeness of Vn{x\z), we obtain the perturbative equation for an{t) as (see Appendix D) 



d 1- Xn\ In 
H Clr, 



dt T J T 



UqJ {2-KYl'^a5ni + Vnan-1 - Vn + la„+i 



1 dz 
2a'dt 



1^ / (n + 2r)! (-1)- 
The center-of-mass position z{t) is determined from the self-consistent condition, 

'^'^ - JZ.dxUix,t)- 

By combining Eqs. ([57j) and ([55]) and denoting the product of the positive integers n, n — 2, . . . 
as nil, we obtain an equation for the bump velocity (see Appendix D), 



dz _2a -^1 + T,n=3,odd V (n-l)!!^" + '^1 , , 

~r — ^ , : — . 

^ UoJi27ry/^a + j:n=o ' 



even V r,n 



Eqs. (|37p and (|39p are the master equations of the perturbative method. For the perturbation 
order specified by n, we include the coefhcients from ao up to o„, and the input components 
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from Iq up to /„, and at least Ii, which is the lowest order of input that drives the dynamics. By 
choosing different orders of perturbation, we can approximate the network dynamics to different 



accuracies. 



The components of the external input in Vnix\z) are calculated to be (see Appendix E) 



In{t) = aUo exp 



8a2 



(40) 



where the noise term rjn{t) satisfies {r]nit)) = and {i]mi't)r]n{t')) = &mn^{t — t'). 
4.1 The Weak Input Limit 

We first compare the perturbative approach with the weak input limit by ignoring all distortion 
modes of the bump; that is, we set a„ = for all n. Due to the readiness of the bump to shift 
its position, the lowest order of the external stimulus to be included is Ii according to Eq. (j39p . 
Thus, in the n = perturbation, Eq. ([39j) becomes 

dz 2a h ^^^^ 



After substituting Ii from Eq. (j40p into Eq. (|4ip . we obtain Eq. (I26p . which is equivalent to the 
network dynamics in the limit of weak inputs. 



4.2 The n = Perturbation 

We further consider the n = perturbation. This amounts to taking into account the effects 
of the change in the bump height on the network dynamics, apart from the positional shift. 
Perturbation up to n = yields 



d 1 - Ao 



dz 
lit 



alio 



(27r)V2aexp 



{zo - zf 



8a2 



r 



aC/o V (2vr)^/2a(zo — z) exp 



+ 2aar]i 



T{Uo^{2TTy/^a + ao) 
With the initial condition oq = at t = 0, the solution of ao from Eq. 



aoit) 



* dt' 



Substituting Eq 



■ exp 



1-An 



it - t') 



alloy (27r)^/2aexp 



(ZQ - zf 

8a2 



reduces to 



+ %(*') 



into Eq. 



we have 



dz 
IE 



a 



(zo - z) exp 



(^0 - zf 
8a2 



+ V2DriA R{t)-\ 



(42) 
(43) 



(44) 



(45) 
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where 



R{t) = 1 + a — exp 

JO T 



-it-t') 



8a2 



+ 



(46) 



and ^o(^) is a stochastic variable given by 

/■* dt' 
?o(i) = / — exp 

JO T 



(l-Ao)(t-t') 



(47) 



i?(t) can be interpreted to be the ratio of the bump height in the presence of the external 
stimulus to its weak input limit. Compared with Eq. (j26p . Eq. (j45p has an extra term R{t)~^. 
R{t)^^ approaches 1 in the limit of weak inputs, but, in general, R{t)~^ < 1, implying that the 
increase in the amplitude of the bump has the effect of slowing down the tracking dynamics. In 
the following two applications, we consider the low noise limit with a = 0. 



4.2.1 Tracking a Moving Stimulus 



Similiar to the derivation of Eq. (I32p . we can determine the dynamical equation for the bump 
lag s after the transient period. From Eq. (j45p . 



dt 



T 



ds as ( s'^ \ N-i 
V exp I ) R{t) , 



(48) 



where R{t) is given by 



'■^ dt' 

R{t) = l + a I — exp 



1-^0/, S{t 



l\2 



'it-t') 



8a2 



(49) 



If the stimulus is not too weak, s reaches its steady state value in a short duration. The integral 
can then be computed precisely, yielding 

„2 



R{t) = 1 + 



a 



exp 



(50) 



1-Ao 

Compared with Eq. (j32p . Eq. (j48p has a correction term, R{t)^^ < 1, representing the effect of 
the bump height distortion. 



4.2.2 Tracking an Abrupt Change in the Stimulus 

The dynamical equation in this case is identical to Eq. (I45p with the initial value z{t = 0) = 0. 
To compute the correction term R{t)~^, however, we need to consider that the bump height is 
initially at a steady value determined by the steady-state response of the network to the stimulus 
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at ^; = 0. This means that Eq. 

aC/oV(2vr)V2a 

ao{t) = 



needs to be replaced by 
(l-Ao)tl 



1-An 



■ exp 



* dt' 



+ aC/oV(27r)V2a / _ exp 

^ Jo T 

and the factor R{t) in Eq. (|l5]) is replaced by 



-it-t') 



8a2 



(51) 



R{t) = 1 + 



a 



1- Ao 



exp 



:i - An)t 



+ a 



* dt' 



T 



■ exp 



^„ {zo-z{t')f 



-{t-t') 



8a2 



(52) 



Indeed, R{t) represents the change in height during the movement of bump. The second and 
third terms respectively contribute to the heights near the initial and final positions (see Fig. 5A). 

Fig. 6 compares the simulation results and the theoretical predictions based on the pertur- 
bation approximation of Eqs. ()45p and ()52p . We see that compared with the weak input limit, 
the result is improved a little bit. 



4.3 The n = h Perturbation 

By increasing the order of perturbation, we get an increasingly better approximation of the 
network dynamics. We observe that the improvement is btter when the order of perturbation n 
increases by 1 to an odd integer. This is because odd distortion modes are antisymmetric, which 
better approximate the bump's shape distortion when the position of the external stimulus is 
biased towards one side of the bump. 

Below, we present the results of the n = 5 perturbation, which correspond to taking into 
account distortions in the position, height, width, skewness and the fourth- and fifth-order 
distortions of the bump. 

By including perturbations up to n = 5, Eqs. ([ST]) and ([39|) become 

d Iq aidz V2 \/6 

^ /^^f/o + ao-^a,U^-#a3 + ^a5, (54) 



dt r V ^ J 2adt 8r 64r 

d l\ I2 V2ai-V3a3dz Vs 

+ 27j = 7 2-a - ^"'^ 

d 3\ I3 \/3a2 - \/4a4 \/5 

+ 4^j = 7 2-a 1-t - 167"^' 

d 7 \ h \/4a3 - \/5a5 dz 

Jt^^r = - 2a dt' 
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2a 1 


[/i + Is + ^/fh + ai) 


1 


r\ 


0(2^) V2a[/o + ao + + y 





(58) 



and as is determined by the condition ai + ^^8/203 + y^lB/Sa^ = 0. Fig. 7 compares the 
theoretical predictions with the simulation results. They agree very well even when the jump 
sizes are large. 



4.4 Shape Distortion 

Besides the reaction time, it is instructive to study the distortion in the bump's shape during 
tracking. As shown in Fig. 8, there is an abrupt change in the height, position, width and 
skewness immediately after the external stimulus is abruptly removed at t = 0. ao{t) approaches 
roughly, as if there were no external stimulus, and there is even a slight overshoot. ai{t) also 
drops abruptly, showing that the peak position lags behind the center of mass, an indication 
that the bump is abruptly distorted. At the same time, distortions in the width and skewness 
increase abruptly, showing that the bump is suddenly pulled by the external stimulus from its 
new position. After the initial abrupt changes, the distorted components undergo rather smooth 
changes on the journey to the new position. Finally, when the bump approaches its destination, 
the height returns to its equilibrium value, while the other components approach 0. 

We also compare the perturbation and simulation results in Fig. 8. The agreement is very 
good. Discrepancies arise from linearizing the interaction term in going from Eq. ([1]) to ([6]). 
Their effects are mainly seen in the moments immediately after the abrupt change in the external 
stimulus. 

Figure 9 shows the distorted bump in the center-of-mass frame during tracking. During the 
process, the lowering of the height, the retarding of the peak position, the broadening of the 
width, and the skewing due to the attraction of the newly positioned stimulus are clearly visible. 



5 Dynamics of a Two-Dimensional CANN 

We can extend our analysis to a two-dimensional(2D) CANN. Consider a neural ensemble 
encoding a 2D continuous stimulus x = (xi,X2), with its components xi G (—00,00) and 
X2 £ (—00, 00). We denote by U{x) the synaptic input to neurons having the preferred stimulus 
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X and r(x) the corresponding firing rate. For a solvable model, we consider the network whose 
dynamics is governed by 

9[/(x,t) 



at 

r(x,t) 



/oo roo 
/ (ix'J(x,x')r(x',t) - C/(x',t), 
-oo J —oo 



' — oo J —oo 



(59) 
(60) 



l + feProoroodx'C/(x',t)2' 

where r is the time constant and p the neural density, and the recurrent interactions are chosen 
to be 



J(x,x') 



27ra2 



exp 



(^x — X 



r\2 



(61) 



where (x — x')^ = {xi — x'^)^ + (x2 — ^2)^ is the Euclidean distance between x and x'. 

It can be checked that when /ext(x, t) = and < < fee = A^p/(327ra^), the network holds 
a continuous family of stationary states given by 

(X-Z)2' 



C/(x|z) 
f(x|z) 



[/q exp 



ro exp 



4a2 
2a2 



(62) 
(63) 



where C/q = [1 + (1 - A;/fcc)^/^]A/(87ra2A:) and ro = [1 + (1 - k/kcy/^]/{A7ra^kp). The parameter 
z = (zi,Z2) is a free variable, indicating the location of the bump. We note the differences in 
J(x,x'), kc, Uq and when compared with their values in the ID case. 

We analyze the stability of the bump states, and consider the network state U{x, t) = 
[/(x|z) +5[/(x, t). After linearizing Eq. ([59|) at {7(x|z), we get 

d 



roo roo 

T—5U{x,t)= / dx'F(x,x'|z)(5C/(x',t) -5?7(x,t) 

Ot J — 00 J —CO 



(64) 



where the interaction kernel is given by 



F(x, x'lz) = —75— exp 



(^x — X 



/\2 



exp 



(x' - 



1 + y/l - k/kc 

27ra2 



exp 



(x - zf 

4a2 



4a2 



exp 



(x' - zf 



4a2 



(65) 



We can construct the right eigenfunction of F{x, x'|z) by using the product of the eigenfunctions 
Un{x\z) in the ID case, i.e.. 



Um,n{x\z) = Um{xi\zi)Un{x2\z2), m,n = 0,1,2, 



(66) 
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The corresponding eigenvalues are given by (see Appendix F) 

Ao,o = Ao (67) 

Am,0 = m / 0, 

Ao,n = A„, n / 0, 
Am,n = AmAn, m^O,n^O, (70) 

where Am are the eigenvalues in the ID case given by Eqs. (jlSp and (I16p . 

The eigenfunction no,o(x|z) corresponds to the amplitude mode of the 2D bump state of the 
network. uifi{x.\z) corresponds to the displacement of the bump position along the coordinate zi 
and tio,i(x|z) the positional shift along Z2- A linear combination of uifi{x\z) and no,i(x|z), i.e., 
cini^o(x|z)+C2tio,i(x, z), with ci and C2 being constants and satisfying c\+c\ = 1, corresponds to 
the positional shift of the bump along the c = (ci, C2) direction. Interestingly, the eigenvalue for 
cini^o(x|z) + C2tio,i(x, z) is 1, indicating that the bump is neutrally stable in the two-dimensional 
space z. 

Analogous to the ID case, we can develop a perturbation method to solve the network 
dynamics, e.g., by using the 2D wave functions, t'm,n(x|z) = Vn{xi\zi)vn{x2\z2) ■, of the quantum 
harmonic oscillators as the basis functions. A more comprehensive treatment of the problem 
involves transforming the basis functions to polar coordinates; it will be presented elsewhere. 
Here, as an illustration, we will present only the result for the weak input limit in the rectilinear 
coordinates, which corrsponds to projecting the network dynamics on a single displacement 
mode. We choose the external input to be 

(X-Z0)2' 



/cxt(x,t) = Qt/oexp 



+ a??(x,t), (71) 



4a2 

where z*^ represents the stimulus location. 

Suppose the bump position is z(t) at time t. Then, the stimulus will pull the bump towards 
the location . In the limit of weak inputs, the network dynamics is dominated by the positional 
shift of the bump along the straight line through z^ . As derived in Appendix F, 



d{z — z^) a ^ 
Jt ~ r 



— [z — z^) exp 



0\2 



z") 



8a2 



+ VZDr?(t), (72) 



where D = a'^ /(Uqt'^tt). This equation describes how the distance between the bump and the 
stimulus varies with time. It has the same form as Eq. (j26p except that z — z^ is replaced by 
z — z^. Thus, the network exhibits the same dynamical properties as observed in the ID case. 
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Fig. 10 compares the theoretical predictions and the simulation results of the reaction time of 
the 2D CANN. As expected, they agree with each other very well with small jump sizes. 



6 Discussion and Conclusions 

The idea that the landscape of the energy function of a CANN has a canyon structure, which fa- 
cilitates the tracking capacity of a neural system has long been recognized (see, e.g., (jAmari 19771 
Hansel and Sompolinsky 1998 Brody, Romo, and Kepecs 2003 IWu and Amari 2005P ). How- 



ever, a detailed and rigorous analysis of this property has been lacking. The main challenge here 
is to solve the network dynamics analytically, elucidating how the shape of the network state 
is affected during the tracking. In this study, by using a simple analytically solvable model, we 
systematically investigate this issue. There are two parts to our main results. First, we show 
that the dynamics of a CANN can be decomposed into different motion modes, corresponding 
to distortions in the height, position, width, skewness and other higher-order distortions in the 
shape of the bump states. We demonstrate that with a CANN, it is the positional shift that 
dominates the network dynamics, with other motion modes having high-order contributions. 
Second, we develop a time-dependent perturbation approach to approximate the network dy- 
namics. Geometrically, this corresponds to projecting the network dynamics onto its dominating 
motion modes. Simulation results confirm that our method works very well. Two interesting 
phenomena in the tracking performance of a CANN are observed, namely, the maximum track- 
able speed and the logarithmic reaction time of the network. Since both properties are associated 
with the unique dynamics of a CANN, they may be tested in practice and serve as clues for 
checking the application of a CANN in neural systems. 

To facilitate an analytical solution and hence to describe the network dynamics clearly, we 
have used a simple firing-rate-based model with the inhibition effect in the form of divisive 
normalization, and the generality of our results may thus be a concern. In particular, we should 
consider whether the mechanism of divisive normalization is plausible, whether tracking behavior 
can be observed in models with different inhibitory mechanisms, and how the perturbation 
method developed in this paper can be used in other types of networks. 

First, we argue that a neural system can have resources to implement divisive normal- 



ization (Grossberg 1988 Heeger 1993 Deneve, Latham, and Pouget 1999). Let us consider a 



network in which all excitatory neurons are connected to a pool of inhibitory neurons. If the 
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time constant of inhibitory neurons is much smaller than that of excitatory neurons, and if these 
inhibitory neurons inhibit the activity of the excitatory neurons in a uniform shunting way, then 
the effect of divisive normalization can be achieved. 

Second, even if there exists uncertainty of the model, the main conclusions of tracking be- 
havior in this work are applicable to general cases. This is because our calculation is based on 
the fact that the dynamics of a CANN is largely dominated by the positional shift of the network 
states, a property coming from the translational invariance of the neuronal interactions, rather 
than from the inhibition mechanism. To see this point clearly, let us consider the dynamics of 
a network having the same form as Eq. ([2]), but the exact form of the inihibition mechansim 
remains unspecified. We consider J[x — x') and r{x) that are properly chosen (e.g., J{x — x') is 
of the Mexcian-hat shape and r{x) a sigmoid function of the inputs), so that the network holds 
a continuous family of stationary states, denoted as U{x\z) and r[x\z), with z representing the 
bump position. This implies that 



oo 



JJ{x\z) = p i dx' J{x — x')f{x'\z). (73) 



By differentiating both sides of the above equation with respect to z, we get 

dU{x\z) f°° , , f°° r „ w /N dr(x'\z) dU(x"\z) 

— pi dx I dx J[x — X )- 



"OO J — oo dU{x"\z) dz 



oo 



dx"F{x,x"\z)^^^^^. (74) 

Thus, the mode of the positional shift of the bump, dU{x\z)/dz, is the eigenfunction of the 
interaction kernel F with the eigenvalue equal to 1. We should observe similar tracking behavior 
of the network, once J{x—x') and r[x) are properly chosen to ensure that all the other eigenvalues 
are less than one. 

Third, the perturbation method developed in this paper can be used to analyze other types 
of networks with tracking behavior. Consider appropriate choices of J{x — x') and r{x) that can 
hold a continuous family of bump states centered at position z (for example, J{x — x') having 
a Mexican-hat shape and r[U) a sigmoid function). Then, the dynamics of the bump can be 
described by Eq. ([6]), with 

F(x, x'\z) = r dx'J{x - ^')^^2^. (75) 
J-oo dU[x'\z) 

Using the basis functions of the quantum harmonic oscillator, we can obtain the matrix elements 
Fmn from Eq. (|13p . While elegant expressions such as those in Eq. (jl4p may not be available in 
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the general case, we can nevertheless compute Fmn numerically. Perturbation dynamics similar 
to Eqs. (I37p and ()39p can be worked out following Appendix D. Indeed, we have considered 
networks with Mexican-hat-shaped interactions. Tracking dynamics with a logarithmic reaction 
time and a maximum trackable speed is observed, and the convergence of the perturbation 
method is equally remarkable. These results will be reported elsewhere. 

The tracking dynamics of a CANN has also been studied by other authors. In particu- 
lar, Zhang proposed a mechanism of using asymmetrical recurrent interactions to drive the 



bump, so that the shape distortion would be minimized (Zhang 1996). Xie et al further pro- 



posed a double ring network model that uses the inputs from the vestibular system in mammals 



to achieve these asymmetrical interactions (Xie, Hahnloser, and Seung 2002). This mechansim 



works well for the head-direction system, but it is unclear if it can be achieved in other neural 
systems. For instance, in the visual and hippocampal systems, it is often assumed that the bump 



movement is directly driven by external inputs (see, e.g., (Samsonovich and McNaughton 1997 



Berry II, Brivanlou, Jordon, and Meister 1999t|Fu, Shen, and Dan 2001[ )), and the distortion of 



the bump is inevitable (indeed the bump distortions in (Berry II, Brivanlou, Jordon, and Meister 1999 



Fu, Shen, and Dan 2001] ) are associated with visual perception) . The contribution of this study 



is that we quantify how the distortion of the bump's shape affects the network tracking per- 
formance, and we obtain the maximum trackable speed of the network. We hope that this 
study will enrich our knowledge about the tracking dynamics of a CANN, which may differ 
in different brain functions. Apart from the tracking capacity, recent experimental observa- 
tions have suggested that CANNs may have other important roles in neural cognitive func- 
tions. For example, the highly structured state space of a CANN may provide a neural ba- 



sis of encoding the categorization relationships of objects (Jastorff, Kourtzi, and Giese 2006 



Graf, Wichmann, Bulthoff, and Scholkopf 2006). It is quite possible that the distance between 



two memory states in the canyon defines the perceptual similarity between the two objects. We 
will explore these issues in the future, and we expect that the theoretical model and approaches 
developed in this work will provide valuable tools. 
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Appendix A: The Hermite Polynomials 

Since the Hermite polynomials are frequently used in the calculations of this paper, here we 
summarize their main properties for the convenience of the reader. 
The nth-order Hermite polynomial is defined as 



The first five Hermite polynomials are given by Hq{x) = 1, Hi{x) = 2x, H2{x) = — 2, 
Hz{x) = — 12x, Hi{x) = 16x^ — 48x'^ + 12. For the analysis in this paper, it is instructive 
to note that the Hermite polynomials satisfy the following properties. 

• Recursion relation: 




(76) 




(77) 




(78) 



(79) 



• Contour integral representation: 




(80) 



with the contour encircling the origin. 



Appendix B: The Interaction Kernel F{x,x'\z) 





(81) 



where we have dropped the dependence on t for brevity, and B = 1 + kp J^^dx'U{x')'^ = 
pAUq/\/^. Interchanging the dummy variables x' and x" in the last term, we obtain Eqs. ([6]) 
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and d?]), with 



F{x, x'\z) 



exp 



1 + - k/h 
lira 



{x — x') 
2^ 



i\2 



exp 



{x' - zf 
4a2 



exp 



{x — z)2 




{x'-zf 


4a2 


exp 


4a2 



(82) 



By substituting the expressions of f7(x), J{x,x') and r(x), we get Eq. ([7]). 
To calculate the matrix elements F^n in Eq. (jl3p . we denote 



(83) 



From Eqs. © and dSD, 



9nix\z) 



-l)"2(^/2a)"- 
a\/7r3/2n!2" 



dx' exp 



[x — x 



/\2 



2a2 



(1 + ^i-k/k,){-inv2ar~^ 



I 



exp 



dx' 

{X — z) 
4a2 



exp 



(x' - Z)2 



2a2 



dx'( A) exp 



(X^ - Z)2 

2^ 



(84) 



The second term vanishes for n > 1. Hence, by evaluating the two integrals explicitly for n = 0, 
we have 



gQ{x\z) = (1 - \Jl- k/kc)vQ{x\z). 



(85) 



Making use of the orthonormality of the basis f„(x|z), we obtain the first line of Eq. ([T 
Next, we consider the case n > 1, focusing on the first term. 
Utilizing the identity (-(i/(iz)"e-(^-^)'/2a2 ^ ^^i^^Y^-{x-z)'^/2a^ ^ ^laYe 



gn{x\z) 



-l)"2(^a) 



d 



dx' exp 



(x - x'y 



aV7r3/2n!2" V dz 
Evaluating the integral explicitly, we find that 

, , , (-l)"2(^/2a)"-5 / dy 
qn(x\z) = , — ' ' 

We are now ready to derive Fmn from gn{x\z): 



2o? 



exp 



ix' - zY 



2a? 



Tx^ 



(x — z)'^ 
4a2 



(86) 



(87) 



(ixi;m(x|z)5(„(x|z). 



From Eqs. P and (|86|) . 

^(-l)'^(^/2a)™-5 (-l)'^(^a)"-^ 



V7rV2m!2"^ \/7rV%!2" 



dx exp 



\x - z) 
4a2 



dx 



exp 



(^x — 



dx 



exp 



(x — z)^ 
4a2 



(89) 
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Using the definition of the Hermite polynomials in Eq. ()76p . this simplfies to 



-lT{y/2a 



m+n— 1 



1 



Vvrm!n!2'"+" (2a) 
Applying integration by parts m times yields 



dxH„ 



F —2 

^ ran — 



V7rm!n!2™+" (2a)" 
Applying Eq. (j78p m times yields, for n > m, 



dxexp 



X — z 
2a 



(x — zY 
2^ 



dx 



exp 



dx) " 



{X — Z) 



X — z 



Fmn — 2 



(\/2a)'"+"-i 1 



Vvrm!n!2™+" (20)*^ a'^(n - m)\ 



POO 


{x-zf 




/ exp 

/ — oo 


2a2 


Hn—m 



2a 



X — z 
2a 



(90) 



(91) 



(92) 



and Fr, 



otherwise. 



Introducing the contour integral representation in Eq. (|80p leads to 



dx exp 



Uc — z 



2a2 



dt 



2"av7rm! 

Interchanging the order of integrating x and we have 



exp 



X — z 



(93) 



F —2 



1-n 




n\ /• at I 
mij 2'Kit'^-'^+^ \ "2 



(94) 



According to the residue theoreom in complex analysis, the contour integral can be reduced to the 
coefficient of in the Laurent expansion of exp(— ^2/2), which is equal to {—l/2)^'^~^^'>^'^/[{n— 
m)/2\\ when n — m is even and otherwise. Hence, we have derived Fmn in the second line of 
Eq. (m. 

Other than the above cases, all other elements Fmn are zero. Hence we have derived Eq. ([1] 
completely. 

As an illustration, we write down the elements of F in the first four rows and columns. 



/ 1 - VI - k/kc 

1 

i 



V 



V2 





8 



...\ 



(95) 
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Since all elements of the matrix F in the lower half of the diagonal line vanish, the eigenvalues 
of F are given by the diagonal elements, namely, Eqs. (jlSp and ()16p . It is straightforward to 
check that the first four right eigenvectors of Fmn are given by 



qo = (1,0,0,...)^ 



(96) 



27 



qi 



(o,i,o,...f, 

^/l72 I - 2^1 -k/k, T 



1 „ /6 „ 



(0,^/-,0,^-,0,...)\ 



(97) 
(98) 

(99) 



where Dq is given in Eq. (j2ip and all the omitted elements are zero. The right eigenfunctions 
of F{x,x'\z), with eigenvalues A^, are given by 



u^{x\z) = ^qnivi{x\z). 
I 

This can be verified by considering the inner product 

/oo roo 
dx'F{x,x'\z)uk{x'\z) = / dx'^v.m{x\z)FmnVnix'\z)^qkivi{x'\z). 
-oo J-oo I 

Making use of the orthonormality of vi{x'\z), we have 

/oo 
dx'F{x, x'\z)u§{x'\z) = Vm{x\z)Fmiqki- 

ml 

Since is an eigenvector of F with eigenvalue A^, 

/oo 
dx'F{x,x'\z)u^{x'\z) = \kYvm{x\z)qkm = Afcuf (x|z). 
-oo ^ 

The first four right eigenfunctions of F are presented in Eqs. (jl7p -(|20p. 
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Appendix C: The Network Dynamics in the Weak Input Limit 

We let U{x,t) = U{x\z{t)) and multiply both sides ofEq. ^ hy vi{x\z{t)), and then we integrate 
over X. Since 

dU{x\z{t)) _ Uo^{2TTy/^adz 



dt 2a dt 

the left-hand side of the resultant equation becomes 



vi{x\z{t)), 



L.H.S. 



TUoJ{27Ty/^adz 



(104) 



(105) 



2a dt 

On the right hand side of Eq.([2|), the first term arises from the external input, consisting of both 
signal and noise components according to Eq. (j25p . The signal component, when projected onto 
the component vi{x\z{t)), leads to 



ii.-ff.S'.signal 



alio 



(27r)i/2 



dx exp 



jx - Zpf 
4a2 



X — z 



exp 



{x - zf 



4a2 



(106) 
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Explicit integration yields 



R.H.S. signal = O-Uq^ (27r)V2fl ^°^^ ^ exp 



(^0 - zf 
8a2 



The noise component leads to 



R.H.S.r 



crr]i{t), 



(107) 



(108) 



where r]i is given in Eq. (|27p and satisfies Eq. (j28p . 

The second and third terms of Eq. ([2]) vanish owing to the fact that U{x\z{t)) is the stationary 
state solution of Eq. ([2]) in the weak input limit, independent of the bump position z{t). 

Combining the above results, we get Eq. (j26p . 

Appendix D: The Perturbation Method 



D.l The Perturbation Equation 

Using Eq. ([8]), we have 



—Vn{x\z) = = 

y^(2^)V2a„!2« 



exp 



(x - zf 
4^ 



X — Z f X — z 



-Hri 



1 rr, '^-Z 



2a? "V^ay V2a "VV2( 



By virtue of Eq. dTTj) and (|78l) . 



Similarly, since U{x\z) = C/oA/(27r)i/2awo(x|z), 



(x|z) = y/ri^Pivn+i{x\z) - ^fri''Vn-i{x\z) 



J_dz 
2a'dl' 



dz 

(109) 



(110) 



dt 



U{x\z) = Uo^i2TT)^/'^aviix\z) 



2a dt' 



(111) 



The neuronal interaction and relaxation terms cancel each other at the steady state. Hence, in 
determining the dynamics, only the linearized terms need to be considered, 

/oo 
dx J{x,x')r{x\z) - U{x\z) ^an^Vm{x\z)Fmn - ^anVn{x\z). (112) 



n m 



Substituting Eqs. (jll0llll2D into Eq. ([2]), we obtain 

^ ^ + ^ 1 ^o\/ (27r)i/2afi + ^ an Vn + Iv^+i - \/nVn-i 



2a 



dz 



(113) 
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Making use of the orthonormality and completeness of Vn{x\z), the coefficients in Eq. (jll3|) 
can be equated term by term. Hence, we get 



dan T 

T— 1 

dt 2a 



By further using Eq. (jl4p . the above equation is re-organized as Eq. ([3 



dz 

-77 = y2 Fnkak - an + In- (114) 

dt ^ 

k 



D.2 The Velocity Equation 

The peak position z{t) of the bump is determined from the self-consistent condition in Eq. (I38p . 
Without loss of generality, we consider 2 = and substitute the expression of Vn{x\z) into 
Eq. (IMl), 

Integrating by parts and using the identity in Eq. ([78|) , we have 

„2 \ 



E 



n=l,odd 



nan 



dxexp ( ) Hn-i 



4a^ 



V2a 



0. 



(116) 



The integral can be performed by using the contour integral representation of the Hermite 
polynomial in Eq. (jSOp . 



E 



n\an 

Vn\¥ J- 



dx (f 



X 



2tx 



n=l,odd 

Interchanging the order of integration yields 



0. 



(117) 
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The contour integral is equal to the coefficient of t" ^ in the series expansion of exp(t^), which 
is [(n - l)/2]!. Hence, 



n=l,odd I ^ ' 

Multiplying both sides of Eq. ([37|l by ^n!!/(n — 1)!!, and summing over odd n, the time deriva- 
tives of a„ disappear, leaving behind 



1 



= 7 E 



nu 



?i=l,odd 



(n- 1)!!^ 



Uoy (27r)V2a + ^ (^^a„_i - Vn + la^+i 



n=l,odd 



1 

2a "dt 



n=l,odd 



'(n-f2r)! (-1) 



(n-1) 



r=0 



n! 



+3r-l».l"'^+2'-- 



(120) 
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Collecting terms in the summations, we have 



;7oV(27r)i/2a+ ^ 

n=0,even 



(n + 1)!! 



(n)!! 



Vn + 1 



'(n- 1)!! 
fn-2)!! 



1 

2a 



1 



7 E 



n!! , \ 1 

— — (/„-a„) + - }^ 



n=l,odd I 



"m\ 



■(m-l)/2 



E 



(-ir 



^ 2m+r--i(^ _ 2r - l)!!r! 



a™.(121) 



m=l,odd 

Note that for odd m, {m — 2r — 1)!! = 2^'^~'^^^^^/'^[{'m — l)/2 — r]!, so that the summation over 
r in the last term is proportional to the binomial expansion of [1 + (_i)](™-i)/2^ which is equal 
to 1 for m = 1, and vanishes otherwise. Therefore, Eq. (|39|) is obtained. 



Appendix E: The Components of the External Input 



We first calculate the component of the signal term in Eq. ([25 

(x - Zof 



In = alio / dx exp 



4a2 



Vn{x\z). 



(122) 



After substituting Eq. ([8]) and completing square in the exponential argument, we have 



alio 



y^(27r)V2an!2'^ 



exp 



(^0 - z)' 
So? 



dx exp 



{x - (ZQ + z)/2f 
2a2 



(123) 



V2a 



Using the identity in Eq. (|79p . denoting [zq + z)/2 as z, and rewriting the integrand in terms of 
the basis functions, we obtain 



In 



y^(27r)V2an!2" 



exp 



(^0 - zf 
8a2 



dxd {2'KYl'^avQ{x\z) 



n 

X E 

A;=0 



k 



V 



{2TTy/^ak\2>'Vk{x\z) , 
I \ y2a 



Zq- Z 



n—k 



(124) 



By the orthonormality of only the = term remains non- vanishing, leading to the signal 
term in Eq. (I40p . 

For the noise component, 



dxr]{x, t)vn{x\z). 



(125) 



Vnit) = 

It is straightforward to check that 

/oo 
dx{r]{x,t))vn{x\z)=0, (126) 
-CXD 

/"CxD /'OO 

{rin{t)r]m{'t)) = dx dx'{rj{x,t)r]{x' ,t'))vm{x\z)vn(x'\z) = 6mnS{t - t'). (127) 



oo J —oo 
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Appendix F: The Dynamics of the 2D CANN 

F.l The Eigenvalues and Eigenfunctions of F(x, x'|z) 
Define 

Mm,n.(x|z) = Um{xi\zi)Unix2\z2). (128) 

We check that u*"'" is the right eigenfunction of F(x, x'|z) and calculate the corresponding 
eigenvalue. 

We distinguish three different cases: 

• m = n = 0. Following the calculation in Eq. (j85p and noting the integrations with respect 
to x'l and X2 can be carried out separately, we have 

/oo poo poo poo 

/ dx'F(x,x'|z)uo,o(x'|z) = / / dx[dx2F{x,x[,x'2\zi,Z2)uo{x[\zi)uo{x'2\z2) 
-00 J —00 J —00 J —00 



= (1- ^1- A://c,)^zo,o(x|z). (129) 

• m 7^ 0,n = (the case for m = 0, n 7^ is the same). An important observation in 
the calculation of Eq. (j84p is that for m > 1, the integral in the second term vanishes. 
Therefore 

/oo poo 
I (ix'F(x,x'|z)iim,o(x'|z) = Amlim,o(x|z). (130) 
-oo J —oo 

• m 7^ 0, n 7^ 0. 

/oo poo 
/ dxF{x,x\z)Um,n{^'\z)=XmKUm,ni^\2')- (131) 
-oo J —oo 

F.2 The network dynamics in the weak input limit 

In the limit of weak inputs, we assume that the shape of the bump is unchanged and only 
the position of the bump varies with time, i.e., U{x,t) ~ f/(x|z(t)). We project the network 
dynamics in Eq. ()59p on a single primary displacement mode, which is the movement of the 
bump along the straight line towards the stimulus location z^. Expressing the left-hand side of 
Eq. (j59p in terms of the basis functions leads to 



a - rC/o(2vr)V2 

T — C/(x|z,t) - 



dzi dz2 



(132) 



at ' ' ' ' 2 

On the right-hand side of Eq. ()59p . we note that the neuronal interaction and relaxation terms 
cancel each other at the steady state, leaving behind the external input. Consider first the signal 
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part, projected onto the components ui^o and uo,i: 



^oo roo 

R.H.Ssi^nni = alio I / dxexp 

-oo J ~oo 



(X — Z 



0^2 



4a2 



[tii,o(x|z) + 'Uo,i(x|z)]. (133) 



Explicit integration yields 
iJ.F.A.signai = -a exp 



For the noise component, 



(z-z0 )2 
8^ 



[(zi - z?)ui,o(x|z) + (Z2 - z^)uo,i(x|t.z)]. (134) 



R.H.S.noise = <Tr?i,o(t)ui,o(x|z) + (T?7o,i (t)tio,i (x|z), 



(135) 



where 



{vMt)mAt')) = {vo,iit)vo,iit')) = 6{t-t'). 



(136) 
(137) 



Using the orthonormality of Um n and expressing the result in vector form, we arrive at Eq. (j72p . 
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Figure 4: (A) The function g{s) and the stable and unstable fixed points of Eq. ([32l) . (B) The 
time dependence of the separation s starting from different initial values, v = 0.025. (C) The 
dependence of the terminal separation s on the speed v. Parameters: N = 200, k = 0.5, a = 0.5, 
T = 1, A = V2Tra, p = N/2tt and a = 0.05. 



38 



A ^ 




Figure 5: (A) The tracking process of the network. (B) The reaction time vs. the jump size 
\zo\. Parameters: 9 = 7r/200 and the rest are the same as Fig. 4. 
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Figure 6: The reaction time vs. the jump size. Parameters: the same as Fig. 5. 
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Figure 7: The reaction time vs. the jump size. The simulation result is compared with the 
n = 5 perturbation approximation. Parameters: the same as Fig. 5. Note that the predictions 
of the n = and n = 1 perturbations are the same because of the self-consistent condition in 
Eq. (|119p . Also, the predicted reaction time given by the n = 1 and n = 2 perturbations are 
almost the same. So are the n = 3 and n = 4 perturbations. 
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Figure 8: The first four a„(t) obtained by the n = 10 perturbation (dashed lines) and by 
projecting U{x,t) onto the corresponding Vn{x\z) during tracking (symbols). Parameters: N = 
200, a = 0.05, a = 0.5, t = 1, k = 0.5, p = iV/27r and A = Vq^tt. 
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Figure 9: Snapshots of the neuronal inputs in the center-of-mass frame of the bump. Symbols: 
simulations. Solid lines: predicted synaptic input with the n = 10 perturbation. Parameters: 
N = 200, a = 0.15, a = 0.5, t = 1, k = 0.5, p = N/2'k and A = V2a7r. 
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Figure 10: (A) The stationary state of the 2D CANN. (B) Illustrating the tracking process of 
the network; (C) The reaction time vs. the jump size. The simulation result is compared with 
the theoretical prediction. Parameters: iV = 40 x 40, = 0.5, a = 0.5, t = 1, A = \/27ra^, 
p = N/{2TTf, a = 0.05 and 9 = V2tt/N. 
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